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DESCRIPTION  AND  PURPOSE 


Potential  application  of  the  stable  laws  has  long  been  hindered  by 
the  unavailability  of  generally  available,  well-documented  algorithms. 

This  paper  removes  this  deficiency  by  presenting  an  algorithm  for  estimation 
of  stable  law  parameters,  with  the  goal  of  facilitating  the  application  of 
stable  laws  in  modeling  and  inference  frameworks.  The  stable  laws  have 
steadily  increased  in  importance  to  the  statistical  community  since  the 
paper  of  Mandelbrot  (1963).  Their  role  as  the  only  laws  possessing 
domains  of  attraction  makes  the  stable  laws  an  appealing  probabilistic 
model,  and  they  are  capable  of  modeling  a  wide  range  of  skewness,  heavy 
tailedness,  and  central  peakedness.  Procedures  for  estimation  of  stable 
law  parameters  have  been  described  by  Mandelbrot  (1963),  DuMouchel  (1971), 
Fama  and  Roll  (1971),  Paulson,  Holcomb,  and  Leitch  (1975),  Koutrouvelis 
(1980,1981),  Feuerverger  and  McDunnough  ( 1981a, 1981b) ,  and  Brockwell  and 
Brown  (1981).  Because  of  the  intractability  of  stable  densities,  attention 
has  centered  in  recent  years  on  Fourier-based  procedures,  using  the 
empirical  characteristic  function.  Such  procedures  should  have  an  adap¬ 
tive  nature  (Paulson,  Holcomb,  and  Leitch,  1975;  Paulson,  Delehanty,  and 
Brothers,  1982;  Paulson  and  Delehanty,  1982). 

IV,  4.V-',r.'-rT 

'  presents  an  iterative  and  adaptive  algorithm  for  joint  estimation 
K 

of  stable  law  parameters,  using  the  empirical  characteristic  function. 

The  algorithm  is  flexible  in  that  either  of  two  procedures  may  be  selected, 
and  subsets  of  the  parameters  may  be  allowed  to  vary  freely,  with  others 
constrained  or  held  constant.  -  The  statistical  rationale  for  the  procedures 
is  described  in  the  companion  paper  by  Paulson  and  Delehanty  (1982).  'The 
algorithm  may  also  be  used  to  provide  informal  estimates  of  the  parameters  ' 
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of  the  stable  law  to  which  a  sample  distribution  is  attracted. 


f<  n-tj  ■: 


f  I 


V  ^  '  ■>  .  i 

J  '  X 

THEORY  AND  NOTATION 

Nondegenerate  stable  random  variables  X  may  be  defined  by  the 
characteristic  function 

$(u)  =  E(expiuX)  =  exp{iuy- |cju  |  a(  1+iB  sgn(u)  x(u>ct))}, 

2 

where  i  =  -1,  0<aS2,|g|sl,  a>0,  and 

X^u,a)  = 


<"  ~ 


(1) 


tan  ,  a/1 


(2) 


^  log|u| ,  0=1. 


Here  a,  the  characteristic  exponent,  is  a  measure  of  heavy  tailedness 

and  central  peakedness,  g  is  a  skewness  measure,  c  is  a  scale  parameter, 

and  n  is  a  location  parameter  unless  (o=l,  0^0),  when  the  function  of 

2 

location  parameter  is  assumed  by  y  +  —  Bologo.  The  only  stable  laws 

IT 

whose  densities  are  expressible  in  closed  form  are  the  Gaussian  (a=2, 

2 

0=0),  the  Cauchy  (a=l,8=0),  and  the  reciprocal  of  a  x  variate  on  one 
degree  of  freedom  (a=J5,0=-l). 

Let  X^,...,Xn  be  a  stable  random  sample.  The  empirical  character¬ 


istic  function  is 


$n(u)  =  n  1  £  exp( iuX. ) 


(3) 


j  =  l 

Let  4»( u )  =  Re  $(u)  +  Im  $(u),  i^nC u )  =  Re  $n(u)  +  Im^tu).  Estimators 
interior  to  the  parameter  space  can  be  viewed  as  zeros  of  the  systems 
Formulation  A 


q  3<|>(u, ) 

l  - L 

3=1 


39 


(ip(Uj)  -  ij;n(Uj))Wj  = 


0 


(4) 


3 


Formulation  B 

1  q  3*(u.)  .. 

I  I  §0“  K3  -  V\))wjwk  =  °*  (5) 

3=1  k=l  J 

for  9  =a,8,a,u.  The  grid  (u^ | j =1,. . . ,q)  is  symmetric  about  zero  but  does 

not  include  the  origin,  and  denotes  the  j  ,k  element  of  the  inverse 

matrix  (K.,  )  1,  where 
3  * 

Kjk =  n  cov<'frn(uj)*  Vuk)) 

=  Re  ♦(u;.-u1<)  +  Im  iKu^+u^)  -  iKu^)  ^(u^).  (6) 

The  weights  (w^ | j=l,. . . ,q}  also  depend  on  the  parameters  a,8,o,u,  and 

are  described  in  the  Numerical  Method  section.  Both  Formulations  A  and 

2 

B  represent  modified,  weighted  x  minimum  procedures,  corresponding 
to  the  respective  objective  functions 


A:  sn  =  I  (*(Uj)  -  *n(u3 ) )2  Wj,  ( 

q  q  -V 

B:  Q  =  I  l  (4»(u. )  -  t  (u. ))  w.  K]K  w.«i(u.  H  (u.  )). 

n  j=l  k=l  ^  "3  3  KKnx 


The  following  points  are  critical  for  practical  application: 

1)  The  shapes  of  i| ;  and  are  highly  dependent  on  location  and 
scale  parameters,  and  so  should  be  standardized; 

2)  The  estimators  are  improved  by  making  the  gridpoints  and  weights 


depend  on  a  and  3; 
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3)  Since  the  procedures  are  adaptive  ( (u. },{w.}  and  {K..  }  depend 

J  J  3  * 

on  unknown  parameters),  algorithms  must  be  iterative; 

4)  Since  a,  2  and  a  are  always  constrained,  each  iteration  in¬ 
volves  solution  of  a  nonlinear  optimization  problem  with 
variable  bound  constraints. 

Our  procedures  may  therefore  be  summarized  as  follows,  where  a  tilde 
indicates  estimators,  their  values,  or  adaptively  standardized 
quantities. 

Begin  with  initial  guesses  for  the  parameters.  At  each  iteration, 
compute  and  save  {u^},  {w^},  possibly  {KJ  },  and  standardized  enpirical 
characteristic  function  values  {^(u.)},  based  on  the  latest  (a, 8).  The 
objective  Sn  or  is  then  minimized  (an  "optimization  subproblem"), 
and  cumulative  location  and  scale  estimates  (£,s)  are  updated.  Iteration 
stops  when  values  of  a  and  u  minimizing  or  are  acceptably  close  to 
unity  and  zero,  respectively. 

Estimators  whose  values  are  not  on  a  bound  are  asymptotically  multi¬ 
variate  Gaussian  distributed.  The  asymptotic  covariance  matrices  £  of 
the  estimators  are  derived  in  Paulson  and  Delehanty  (1982).  The  basic 
formula  is 


T.  =  Hi1  V.  hT1,  i  =  A,B.  (C 

There  are  two  particularly  appealing  ways  to  approximate  £.  In 
"approximation  (i)",  expectations  are  approximated  from  the  data:  H  is 
computed  by  differencing  the  objective  at  the  final  optimum,  and 

VA  =  4DTKD,  (3 

yB  =  4dtk-1k(w)k"1d. 


(11) 
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attracted  to  a  stable  law  with  parameters  (a,8,a,y),  the  sequence  of 
normalized  estimators  a^/k^01^,  u^/k),  for  reasonable  values  of 

k,  should  approach  (o,8*c.u).  In  particular,  a  rapid  rise  in  {a^}  may 
indicate  that  a  stable  model  is  not  appropriate,  a  possible  alternative 
being  a  mixture  of  finite  variance  components  with  differing  scale  para¬ 
meters  . 

The  k-sum  procedure  can  thus  be  used  in  a  sensitivity  analysis,  to 
examine  how  well  the  data  support  the  stability  assumption.  Other  possible 
tools  for  sensitivity  analysis  are  varying  the  mechanism  (to  be  described 
below)  underlying  the  weights  {w^},  and  comparing  approximations  (i) 
and  (ii)  of  the  estimated  asymptotic  covariance  matrix,  provided  n  is 
large  enough  for  approximation  (i)  to  be  accurate. 


NUMERICAL  METHOD 

The  main  computational  task  required  is  solution  of  bound- constrained 
nonlinear  optimization  problems.  Although  Formulations  A  and  B  lead  to 
nonlinear  least  squares  problems,  current  algorithms  for  nonlinear  least 
squares  do  not  allow  constraints  (Hiebert,  1981).  Numerical  Algorithms 
Group  (NAG)  subroutine  E04KBF  (NAG,  1981)  is  used  for  optimization. 

E04KBF  is  a  quasi-Newton  procedure,  requiring  an  objective  function  and 
analytical  first  partial  derivatives.  It  is  substantially  faster  than 
the  gradient  projection  routine  used  by  Paulson,  Holcomb  and  Leitch  (1975), 
although  the  latter  is  very  reliable.  The  other  complicated  numerical 
procedure  required  is  inversion  of  a  positive  definite  symmetric  matrix 
(K,  H  or  H),  for  which  NAG  subroutine  F01ABF  is  used.  Various  NAG  utility 
procedures  are  also  used,  see  Auxiliary  Algorithms.  The  use  of  NAG 
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procedures  inhibits  transportability  in  that  the  algorithm,  as  presented, 
is  only  usable  at  installations  having  the  NAG  Library.  However,  listings 
of  rapid,  high-quality  algorithms  for  constrained  optimization  have  not 
appeared  in  the  literature  (see  Chambers,  1977,  pp.  159-160;  the  situation 
described  there  has  not  improved).  Given  that  E04KBF  is  used,  reliance 
on  additional  NAG  Library  procedures  is  expedient. 

We  require  a  minimum  of  q=20  gridpoints  {u^},  and  prefer  q=20  or 
40,  since  they  are  reasonable  values  in  practice,  and  have  been  tested 
extensively.  Only  the  positive  gridpoints  are  explicitly  required,  due 
to  symmetry  of  the  grid  and  the  Hermitian  property  of  characteristic 
functions.  They  are  computed  as  follows:  An  endpoint  U  is  chosen  as 
3,  c  2  1.8;  3.3,  1.8  >  a  >  1.7;  3.6,  1.7  >  a  >  1;  5,  a=l;  4,  1  >  a  a  .9; 

5,  .9  >  a  i  .8;  7,  .8  >  a  2  .6;  10,  a<  .6.  An  inner  number  I  of  points 
is  selected  close  to  the  origin:  1=2  if  q<30  and  3  if  qa30,  1  being  sub¬ 
tracted  if  as. 5.  The  inner  I  points  are  spaced  as  follows:  if  a>l, 

*5  the  "a-optimal"  values  of  Feuerverger  and  McDunnough  (1981b)  for  the 
nearest  (larger)  a  are  used;  if  aSl,  the  first  I  points  giving  q/2  equal 

Jm  3  *2 

increments  of  log  (u  +  a  )  between  0  and  U  are  multiplied  by  %  a 

A 

(a  =  max(a,  .3)).  The  remaining  points  are  logarithmically  spaced  out 

-  „  3 

to  U:  if  a>l,  the  function  log  (1  +  u/2)  is  used,  and  if  asl,  log(u  +  a  ) 

is  used.  This  rather  complicated  ad_  hoc  scheme  was  developed  through 
graphical  inspection  of  d>(  u )  and  iJi^Cu),  comparisons  of  asymptotic  effic¬ 
iencies,  and  parameter  estimation  for  real  and  simulated  data.  No 
claims  of  optimality  are  made,  but  the  scheme  provides  high  efficiencies 
if  efficiency  is  preferred,  or  good  matches  between  and  i  if  curve 
fitting  is  preferred.  The  point  of  stratified  and  logarithmic  spacing  is 


a 


to  emphasize  u  values  near  the  origin.  Details  when  a£l  reflect  the 
fact  that  t|>(u)  has  a  sharp  cusp  near  the  origin,  but  decays  slowly. 
The  stepwise  nature  of  the  scheme  is  not  deemed  a  serious  drawback. 

The  weights  {w_.}  are  computed  as  follows: 

Under  Formulation  A, 


w . 
] 


U(u^) 


i  2X 


vyV 


exp(-2lju..  |a) 
KT(Uj’Uj} 


(17) 


and  under  Formulation  B, 

Wj  =  |  ( u^  )  ]  X  =  exp(-X  |  u^  |°) ,  (18) 

where  X  and  t  are  supplied  by  the  user,  0  s  t  s  1, 

K^(  u,u)  =  1  +  r(Im  $(2u)  -  ip  2(  u) ) ,  (19) 

and  X  is  recommended  nonnegative.  Rationale  for  these  weights,  and 
some  corresponding  asymptotic  efficiencies,  are  in  Paulson  and  Delehanty 
(1982).  We  recommend  x=l  under  Formulation  A.  Under  Formulation  B, 
it  is  convenient  to  let  x£0  represent  a  fraction  of  the  average  diagonal 
element  by  which  to  inflate  K,  giving  a  matrix  we  shall  call  A.  We 
have  only  found  this  inflation  necessary  if  a  is  very  close  to  two, 
when  x=0.01  suffices. 

To  use  the  quantity  X  as  a  tool  for  sensitivity  analysis,  we  inter¬ 
pret  it  as  a  damping  factor,  lessening  the  effects  of  noise  in  \|>  (u)  for 
larger  |u|.  If  the  data  are  truly  stable  and  the  sample  size  is  fairly 
large  (say  150  or  more),  estimates  should  change  little  as  X  varies, 
say,  from  0  to  1.  Large  discrepancies  in  the  estimates  for  different 
values  of  X  indicate  problems  with  the  data  or  the  stable  assumption  or 
both.  It  may  not  be  easy  to  isolate  the  difficulty  but  further  study  is 


9 


definitely  required. 

For-  the  k-sum  procedure,  k>l,  the  situation  regarding  gridpoints 

and  weights  changes.  Tests  so  far  indicate  that  when  ct>l,  Formulation  A, 

with  gridpoints  equispaced  from  0  to  U,  gives  better  results  than 

"efficient  configurations"  used  for  k=l.  Reasons  for  this  are  unclear. 

*•  *p  k 

A  possible  explanation  is  that  when  a>l  and  k>l,  ^(u)  is  so  smooth  that 

estimation  is  practically  equivalent  to  deterministic  curve  fitting,  and 

implicit  or  explicit  emphasis  on  gridpoints  near  the  origin  neglects 

important  curvature  for  large  |u|.  Accordingly,  when  k>l  and  ot>l,  we 

~k 

equispace  gridpoints  and  set  all  weights  to  1.  When  aSl,  ^(u)  has  a 

sharp  cusp  near  the  origin  and  remains  a  jagged  curve  as  k  increases, 

due  to  the  presence  of  very  large  observations.  In  this  situation,  we 

set  all  weights  to  1  and  use  basically  the  same  gridpoint  scheme  as  when 

*2 

k=l,  omitting  only  multiplication  of  the  inner  points  by  a  /4.  In  either 
case.  Formulation  A  is  recommended. 

An  important  question  is  how  large  k  may  be  taken.  Equation  (16) 
suggests  that  we  cannot  expect  to  take  k  arbitrarily  large.  There  seems 
to  be  a  tendency  for  a  to  increase  and  B  to  drift  if  k  is  too  large, 
though  this  may  be  partially  due  to  suboptimal  gridpoints  or  weighting. 

It  appears  that  when  n  is  large,  say  500  or  more,  k  may  safely  be  taken 
up  to  20.  Care  is  required  for  smaller  n,  and  when  a  is  small  or  very 
near  two. 

Implicit  standardization  is  carried  out  as  follows.  Let  k  be  a 
positive  integer,  and  (£,s)  cumulative  location  and  scale  estimates. 


Then 
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<i(ys  pjk(cos  Yjk +  sin  Yjk}’  (2o) 

where 

Pjk  =  Un(^/s)|k  (21) 

and 

y..  =  k  arg*  (u./5)  -  2u. /s.  (22) 

j  K  lx  J  J 

No  problems  of  principal  values  arise,  and  complex  arithmetic  is  not 
used-  The  FORTRAN  mathematical  library  function  ATAN2  computes  argu¬ 
ments. 

The  estimator  a  may  be  bounded  in  (closed)  subintervals  of 
[6,1-e],  [1,1],  or  Cite, 2]  unless  8  is  fixed  at  0,  when  [6,2]  is  possible 
(5  and  e  are  small  positive  numbers),  while  S  may  be  bounded  in  sub¬ 
intervals  of  [-1,1].  Estimators  c  and  u  may  be  constrained  arbitrarily 
in[6,“)  and  (-«,<»),  respectively,  unless  a  is  fixed  at  1  and  8  is  not 
fixed  at  0,  when  u  cannot  be  constrained.  Bounds  on  a  and  u  are  internally 
set  for  use  in  subproblems.  These  bounds  must  be  wide  enough  to  allow 
the  "true  values"  to  be  found,  but  narrow  enough  to  deter  straying  into 
undesirable  regions,  particularly  a-**,  The  ad  hoc  bounds  of 

[-5,5]  for  u  and  [0.2,5]  for  a  work  well  in  practice.  If  o  or  p  are 
initially  constrained,  their  internal  bounds  are  adaptively  modified,  see 
the  Algorithm  for  description. 

Initial  guesses  for  the  parameters  are  required.  We  do  not  find 
their  specification  particularly  important,  provided  a  is  on  the  correct 
side  of  1  in  the  nonsymmetric  case.  We  have  used  the  median  and  semi- 
interquartile  range  as  guesses  for  u  and  o,  and  averages  of  upper  and 
lower  bounds  for  a  and  8.  If  a  is  anticipated  less  than  1.2,  say,  it  is 
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worthwhile  to  put  more  effort  into  initial  guesses,  since  fewer  iter¬ 
ations  will  be  required  (the  semi- interquartile  range  will  overestimate 

a,  and  if  a  is  near  but  different  from  1,  the  median  is  nearer 
2 

u  -  —  6  a  log  a  than  u). 

Convergence  is  judged  by  a  tolerance  on  subproblem  solutions, 
max(|a(m)-l|,|y(m)|),  (or  max  (jc.(m)  -  c(ra_1) | .  |s(m)  -  8(®"1)|)  if 
5  and  u  are  fixed),  with  a  maximum  allowable  number  of  iterations. 
Attainable  tolerances  depend  on  n,  but  more  strongly  on  the  underlying 
parameters.  If  a  is  near  two,  ij»n  is  very  smooth  and  stringent  tolerances 
such  as  10  6  may  be  attained.  If  asi.2,  has  many  small  oscillations 
due  to  large  observations,  and,  especially  for  smaller  samples,  it  may 
be  preferable  to  terminate  after  a  fixed  number  of  iterations.  Good 
estimates  are  generally  obtained  within  five  iterations,  fewer  if  initial 
guesses  are  good;  if  stringent  tolerances  are  required,  or  for  difficult 
problems  (skewed  distributions  with  0.9SaSl.l)  more  may  be  required. 
Convergence  is  typically  slower  under  Formulation  B,  since  the  weighting 
mechanism  is  more  conplicated. 

Approximation  of  asymptotic  covariance  matrices  requires  little 
description.  We  note  that  for  approximation  (i)  and  the  q  values  we  use, 
it  is  faster  to  define  a  vector 


<5 j  =  (ijku^  -  cos  u^ 


sin  u^  , 


♦(u  )  -  cos 

q 


u  X, 

q  j 


sin  u^ ) , 
(23) 


and  cumulate 

V 


l  <2TV 

j-1  1 


(DT6^)1 


(2U) 


under  Formulation  A,  or 
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V  =  4n_1  l  (DTA'16.)(DTA“15.)T  (25) 

j=l  3  3 

under  Formulation  B,  than  to  cumulate  K.  The  matrix  D  is  computed  by 
the  function/gradient  subroutine.  E04KBF  returns  an  approximate  Hessian, 
which  could  conceivably  be  used  for  H  in  approximation  (i).  Rather 
often,  however,  E04KBF  will  terminate  with  its  failure  indicator  set  to 
3  and  the  Hessian  set  to  the  identity  matrix,  even  though  the  optimum 
may  be  reliable.  It  is  therefore  simpler  to  compute  H  by  differencing. 

The  following  procedures  is  used:  Set  an  initial  Hessian  to  0,  and  the 
steplength  to  10  '3.  Successively  divide  the  steplength  by  /lo’  and  approx¬ 
imate  the  Hessian  by  differencing;  three-point  differencing  for  the 
diagonal,  and  four-point  for  off-diagonal  elements.  Compare  elements  of 
successive  approximations  by  maximum  relative  or  absolute  differences, 
according  as  the  element  of  the  latest  approximation  exceeds  1  in  absolute 
value  or  not.  A  tolerance  of  10-  is  used  for  this  convergence  criterion. 

If  convergence  has  not  occurred  with  a  steplength  of  10  5,  the  result 

.  -4  . 

with  steplength  10  is  used. 

Approximation  (i)  of  the  asymptotic  covariance  matrix  is  rather 
expensive  to  compute.  It  should  not  be  computed  for  smaller  sample  sizes, 
as  it  implicitly  involves  estimation  of  4q(Hq+l)  covariances. 

Following  is  an  informal  description,  in  Algorithm  form,  of  the 
basic  routine  STABLE.  Approximate  asymptotic  covariance  matrices  may 
also  be  computed,  but  this  presents  no  logical  difficulties,  so  is  omitted. 
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Algorithm 

Produces  estimates  (a,0,o,y)  for  k-sums  (k^l),  based  on  a  sample 

( x.  i  •  •  ■ » x  )  • 

1  n 

Input  parameters:  k,  n,  {x^},  q,  X,  t,  Formulation  (A  or  B), 
convergence  tolerance  e,  maximum  number  M  of  iterations,  and  flags 
whether  a  and  u  are  constrained. 

Input/output  parameters:  (a,fi,a,y)  are  initial  guesses  on  entry 
and  estimates  on  exit,  (a,  ,B,  ,uT )  end  (a  ,6  ,a  ,y  )  are  lower  bounds. 

LLLiL  U  U  U  U 

The  (a,y)  bounds  are  changed,  but  restored  on  exit.  In  the  special  case  where 

-  -  2 
a  is  fixed  at  1,  y,  on  entry,  is  the  initial  guess  for  location  u  +  —  Bologo. 

Auxiliary  quantities  t  and  s  are  cumulative  location  and  scale 

estimates.  Entry  values  of  (a^jO^y^y^  dre  St°rec*  ^bl’b2,l53,^4^ '  0n 

entry  and  exit,  (a,y,crT ,uT ,a  , ,u  )  are  normalized. 

u  L  U  U 

Si  [Initialize.] 

**  -  1/a.  - 

Sett*-  ku,  s  +  k  a,  nr-O. 

Save  (b1>i>2»b3 >b4)  ^0l,Ou’UL,uu^ * 

if_  ii  is  unconstrained  then  set  y^  *■  -  5,  y  5; 

else  if  y  is  fixed  then  set  y.  yt  -*-0; 

else  set  y^  kyL,  yy  ky^. 

if  3  is  unconstrained  then  set  aL  0.2,  ou  5; 

else  if  a  is  fixed  then  set  a.  a  •*-  1; 

■  u 

else  set  aL  «-  k1//a  <?L,  au  k1^  au> 
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[Looping  point  for  iteration;  save  adaptive  quantities  for  sub- 
problem.  j 

Increment  m  ■«-  m  +  1. 

0  -(m-1)  -  3(111-1)  3 
Save  a  *■  a,  0  *■  6. 

if  a  is  constrained  but  not  fixed  then  set  or  *■  max  (0.2,b./s), 

a  *■  min  (5,b„/s). 
u  2 

if  y  is  constrained  but  not  fixed  then  set  y  *■  max  (-5,(b_-£)/s), 
— 1  "  Li  0 


yu  +•  min  ( 5 ,(b4-£)/s) . 


Set  a«-l,  y-*-0. 


Compute  and  save  positive  gridpoints  {u^  j  j  =  q/2+1, . . . ,q}, 

weights  (Wj | j  =  l, . . . ,q),  and  standardized  empirical  characteristic 
^  k 

function  values  {^(u^)  |  j=l,...,q). 
if  Formulation  B  then  compute  and  invert  A. 


[Subproblem. ] 

Sclve  the  optimization  problem 


min  1  (<l>(  u.)  -  *k(u.))2  w. 

j=l  1  n  3  3 


(Formulation  A) 


!  \  wi(|<'(ui)  -  $£(1^))  Ai3( 


i=l  j=l 


4<(Uj)  -  in(u^))w^  (Formulation  B), 


yielding  new  ( a, 3, a,y). 


S4  [Update  and  test  convergence.] 


Set  t  *■  t  +  s  5. 


if  a=l  then  set  2  s  S  a  log 
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Set  s  +  s  o. 

if  a  and  u  are  fixed  then  error  max  ( | ,  1 6-S^ m_ 1 ^ | ) ; 

else  error  *■  max  ( 1 5- 1 1 ,  j  w  | ) . 
if  error  a  e  and  m<M  then  go  to  S2_. 

S5  [Final  estimates.] 

Set  (cL’°u’pL,liu^ 

if  a=l  then  set  Z  Z  -  —  §  s  log  s. 

"  7T 

Set  u-*-Z/k,  5+^/k1^. 
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STRUCTURE 


SUBROUTINE  STABLE  (X,N, MODE, KSUM,XLAM,TAU,NPTS,TOL,MAXIT, XL, XB,XH, NPAR, 
ISCLB  D , LOCBND , ICOV, VCV1 , VCV2 , WORK, LWORK, IWORK , L IWORK , IFAULT ) 


Formal  parameters 


X 

Real  array  (N) 

input : 

N 

Integer 

input : 

MODE 

Integer 

input : 

KSUM 

Integer 

input : 

XLAM 

Real 

input : 

TAU 

Real 

input : 

NPTS 

Integer 

input : 

TOL 

Real 

input : 

MAXIT 

Integer 

input : 

XL 

Real  array  (NPAR) 

input : 

output: 

XH 

Real  array  (NPAR) 

input : 

output : 

NPAR 

Integer 

input : 

ISCLBD 

Integer 

input : 

sample 
sample  size 

formulation;  if  zero,  then 
Formulation  B  is  used,  else 
Formulation  A 

convolution  power  k 

X 

T 

q 

convergence  tolerance 

maximum  allowable  number  of  iterations 

lower  bounds  for  parameters  ; 

the  third  and  fourth  elements  change 
during  execution 

input  values  are  restored 

upper  bounds  for  parameters ;  the 
third  and  fourth  elements  change 
during  execution 

input  values  are  restored 

number  of  parameters  (4) 

flag  if  a  is  constrained:  if 
negative,  o  is  fixed  at  XB(3);  if 
zero,  o  is  free  to  vary  and  initial 
values  of  XL(3)  and  XH(3)  are  irrele¬ 
vant;  if  positive,  a  is  constrained 
in  C  XL(  3)  ,XH(  3)  ] 
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LOCBND  Integer  input: 


ICOV  Integer  input : 


VCV1  Real  array (NPAR,NPAR)  output: 


VCV2 

Real  array(NPAR,NPAR) 

output: 

WORK 

Real  array  (LWORK) 

workspace 

output : 

LWORK 

Integer 

input : 

IWORK 

Integer  array( LIWORK) 

workspace 

output: 

LIWORK 

Integer 

input : 

IVNIT 

Integer 

input : 

flag  if  u  is  constrained:  if 
negative,  u  is  fixed  at  XB(4);  if 
zero,  u  is  free  to  vary  and  initial 
values  of  XL(4)  and  XH(4)  are  irrele¬ 
vant;  if  positive.,  u  is  constrained 
in  [XL( 4) ,XH( 4) ] 

flag  for  computation  of  covariance 
matrices:  if  negative,  neither 
approximation  (i)  nor  (ii)  is  computed; 
if  zero,  both  approximations  are  computed; 
if  positive,  only  approximation  (ii)  is 
computed 

covariance  matrix  approximation  (i) 
if  requested;  the  strict  lower 
triangle  contains  correlations,  the 
upper  triangle  contains  covariances 
(times  n);  if  a  parameter  is  on  a 
bound,  the  corresponding  elements  are 
zero 

covariance  matrix  approximation  (ii) 
if  requested;  the  strict  lower  triangle 
contains  correlations,  the  upper  triangle 
covariances  (times  n);  if  a  parameter  is 
on  a  bound,  the  corresponding  elements 
are  zero 


some  elements  may  be  of  interest  on 
output  (see  Restrictions) 


some  elements  may  be  of  interest  on 
output  (Tee  Restrictions) 


if  positive,  unit  number  for  output 
(see  Additional  Comments);  if  zero 
or  negative,  no  output  is  produced 


I FAULT  Integer 


output: 


failure  indicator 


18 


Failure  indicators 

IFAULT  =  0  indicates  success.  Nonzero  values  of  IFAULT  are  due 
to  two  types  of  errors.  The  first  type  is  input  errors,  detected  in  STABLE 
IFAULT  will  be 

1  if  MAXITSO ; 

2  if  N<50  (see  Restrictions); 

3  if  KSUM  S  0; 

4  if  t<0  or  t  >1  and  MODE^O; 

5  if  NPTS <20  or  mod(NPTS,2)*0; 

6  if  T0LS0; 

7  if  NPAR/4; 

8  if  insufficient  workspace  was  allotted  (see  Restrictions); 

9  if  improper  bounds  were  supplied.  The  following  conditions 
cause  this  failure: 

XL(i)>XB(i)  or  XL(i)>XH(i)  or  XB(i)>XH(i),  i=l,2 
XL(1)S0  or  XH( 1) >2 
XL( 1) <-l  or  XH( 1 )  >1 

(XL(2)*0  or  XH( 2)^0)  and  (XL(1)<1  and  XH(1)*1 
or  XL(1)S1  and  XH(1)>1) 

XB( 3)^0 

ISCLBD^O  and  (XL(3)>XB(3)  or  XL(3)>XH(3)  or  XB(3)>XH(3)) 

ISCLBD<0  and  XL(3)*XH(3) 

LOCBND^O  and  (XL(4)>XB(4)  or  XL(4)>XH(4)  or  XB(4)>XH(4)) 

L0CBND<0  and  XL(4)?*XH(4) 

LOCBND^O  and  XL(  1)  =  XH(  1)=1  and  (XL(2)*0  or  XH(2)*0). 

On  input  errors,  STABLE  terminates  immediately,  without  performing  any 
computations.  The  second  type  of  error  occurs  after  some  computation. 
IFAULT  will  be 
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if  A  was  found  numerically  non  positive  definite; 


11  if  A  was  found  ill-conditioned; 

12  if  too  many  function  evaluations  were  required  during  solution 
of  a  subproblem; 

13  if  iteration  converged,  but 

the  most  recent  E04KBF  fault  indicator  was  3  and 
internal  checks  were  not  met.  These  checks  are 
(i)  || G j| 2  <  10  *  X02AAF( DUMMY),  and 
(ii)  K  <  l/|lG||  ,  as  recommended  by  E04KBF 
documentation,  where  ||  G  ||is  the  norm  of  the  projected 
gradient  and  K  the  estimated  condition  number  of  the 
projected  Hessian  matrix; 
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14  if  there  were  repeated  problems  with  overflow  in  the 
Cholesky  factors  of  the  projected  Hessian; 

15  if  iteration  converged,  but 

the  most  recent  E04KBF  fault  indicator  was  5  and 
internal  checks  were  not  met; 

16  if  convergence  did  not  occur  in  MAXIT  iterations ; 

17  if  convergence  did  not  occur  in  MAXIT  iterations,  the 
most  recent  E04KBF  fault  indicator  was  3,  and  internal 
checks  were  not  met; 

18  if  convergence  did  not  occur  in  MAXIT  iterations,  the 
most  recent  E04KBF  fault  indicator  was  5,  and  internal 
checks  were  not  met. 

Conditions  IFAULT=10  and  11  are  detected  in  SETECF  (they  are  caused  by 
t being  too  small  under  Formulation  B),  the  remainder  in  STABLE. 

IWORK( 2 )  and  IW0RK(3)  (see  Restrictions)  are  failure  indicators 
for  asymptotic  covariance  matrix  versions  (i)  and  (ii)  respectively.  Zero 
indicates  success,  1  that  H  was  non  positive  definite,  and  2  that  H  was 
ill-conditioned,  the  failures  detected  in  SETVCV.  If  IFAULT=1-12  or 
14,  covariance  matrices  are  not  computed,  and  their  fault  indicators  are 
set  to  the  corresponding  value  of  IFAULT. 


Auxiliary  algorithms 

The  user  has  only  to  call  STABLE.  Auxiliary  procedures  fall  into 
two  groups:  those  supplied  here,  and  NAG  Library  procedures.  The 
following  subroutines  are  supplied: 

SUBROUTINE  GRI0WT(PAR,NPAR,XLAM,TAU,PTS,NPTS2 ,WT,NPTS ,M0DE ,KSUM) :  computes 
gridpoints  and  weights; 

SUBROUTINE  CHARFN(U,PAR,NPAR,RE,XIM):  computes  real  and  imaginary  parts 
of  standard  stable  characteristic  function  $(u); 

SUBROUTINE  FUNCT( IFLAG,N,XC ,FC ,GC,IW,LIW,W,LW) :  objective  function/ gradient 
evaluation; 
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SUBROUTINE  SETECF(X,N, PAR, NPAR, MODE, TAU, SIGMA, XMU,KSUM,IA,NPTS2,NPTS, 
PTS,ECF,A,AINV,WORK, IFAULT) :  computes  standardized  empirical  character¬ 
istic  function  values  ^(u),  computes  and  inverts  A  under  Formulation  B; 

SUBROUTINE  MONIT(N,XC,FC,GC,ISTATE,GPJNRM,COND,POSDEF,NITER,NF,IW,LIW, 
W,LW):  monitors  the  progress  of  E04KBF; 

SUBROUTINE  VARIAB(  IC0V,X,N , PAR, NPAR, MODE .SIGMA ,XMU,ISUB,NVAR,PTS ,NPTS2 , 
WT,ECF,NPTS,DERIV,W0RK,H0LD,A,IA,AINV,VCV1,VCV2,H,NVAR1,V,IW,LIW,W,LW, 
IFAIL1,IFAIL2) :  computes  approximate  asymptotic  covariance  matrices; 

SUBROUTINE  VMATRX  (X,N, MODE, XMU, SIGMA, PTS,NPTS2,WT,ECF,W0RK,NPTS,DERIV, 
V,HOLD,NVAR) :  computes  7  for  version  (i)  of  asymptotic  covariance  matrix; 

SUBROUTINE  DAPR0D(FAC1,IFAC1,NPTS,FAC2,W0RK,NVAR):  auxiliary  matrix 
multiplication  for  VARIAB : 

SUBROUTINE  HVPR0D( FAC1, IFAC1,NVAR,FAC2 ,NPTS , VH.IVH) :  auxiliary  matrix 
multiplication  for  VARIAB; 

SUBROUTINE  SETVCV( ISUB,NVAR,H, NVARl.V, WORK, VCV, NPAR, SIGMA, IFAULT) : 
auxiliary  routine  for  VARIAB; 

SUBROUTINE  HESDIF( PAR, NPAR, ISUB ,H,SAVE1,SAVE2 ,NVAR, IW,LIW,W,LW) :  computes 
an  approximate  Hessian  by  differencing  for  version  (i)  of  asymptotic  co- 
variance  matrix. 

The  following  NAG  Library  procedures  are  used: 

REAL  FUNCTION  X02AAF( DUMMY) :  returns  the  smallest  positive  e  such  that 
1.0  +  e  >  1.0; 

SUBROUTINE  E0 4KB F( N , FUN CT , MON I T , IPRINT , LQ CS CH , INT YPE , M INL IN , MAXCAL , ETA , 
XTOL,STE?MX,FEST,IBOUND,BL,BU,X,HESL,LH,HESD,ISTATE,F,G,IW,LIW,W,LW,IFAIL) 
solves  optimization  problems.  Control  parameters  are  set  as  follows: 


IPRINT 

= 

0 

LOCSCH 

S 

.TRUE. 

INTYPE 

= 

3  for  subproblems  after  the  first 

if 

not  fixed  are  not  on  bounds ,  else 

0 

MINLIN 

s 

NAG  Library  routine  E04LBS 

MAXCAL 

= 

400 

ETA 

= 

0.9 

XTOL 

= 

10 ,0/X02AAF( DUMMY)  explicitly,  so 

it 

STEPMX 

= 

0.25 

FEST 

= 

0.0 

IBOUND 

= 

0  ; 

SUBROUTINE  FQ1ABF( A,IA,N ,B , IB ,Z, IFAIL) :  inverts  the  positive  definite 
symmetric  matrix  A; 
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SUBROUTINE  F01CAF(A,M,N,IFAIL) :  sets  matrix  A  to  zero; 

SUBROUTINE  F01CMF(A,LA,B ,LB,M,N) :  copies  elements  of  matrix  A  into 
matrix  3; 

SUBROUTINE  F01CKF(A,B,C,N,IP,M,Z,IZ,IOPT,IFAIL) :  matrix  multiplication 
A=BC,  where  B  or  C  may  be  overwritten. 


RESTRICTIONS 

We  require  the  sample  size  N  at  least  50,  since  for  smaller  samples 
ij)n(u)  is  not  generally  sufficiently  smooth  to  allow  accurate  estimation. 
Since  a  and  B  are  bounded  in  the  narrow  ranges  (0,2]  and  [-1,1]  and  have 
standard  errors  decreasing  as  N  ,  it  is  preferable  to  have  N£100.  For 
N  less  than  150,  say,  relatively  large  values  of  X  may  be  preferred,  to 
damp  out  noise  in  iji^Cu).  We  further  require  NPTS£20. 

Extended  work  vectors  WORK  and  IWORK  are  required,  in  order  to 
communicate  information  to  FUNCT  and  MONIT  without  using  COMMON  blocks. 

To  aid  readers  who  may  wish  to  adapt  the  algorithm  to  installations  not 
having  the  NAG  Library,  we  describe  the  use  of  these  work  vectors. 

The  required  length  of  WORK  is  10  +  11*NPAR  +  NPAR*(NPAR-l)/2  + 

( 3+NPAR)*NPTS  +  NPTS  +  NPTS/2  if  MODE*0,  with  an  additional  NPTS*(2*NPTS+1) 
required  if  M0DE=0.  Some  sample  lengths  are 

MODE  N?TS=20  NPTS=4Q 

0  10  30  3600 

nonzero  210  360 


The  subvector  W  is  passed  to  E04KBF, FUNCT,  and  MONIT. 


WORK  starting  point  W  starting  point  Elements  Used  for 


\ 
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The  required  length  of  IWORK  is  7  +  NPAR.  The  subvector  IW  is 
passed  to  E04KBF,  FUNCT,  and  MONIT  . 


IWORK  starting  point  IW  starting  point  Elements  Use  for 


1 

2 

(other  addresses 
internally 
computed) 


1 

3 

4 

5 

6 


1 

NPAR 


2 

1 

1 

1 

1 


Iteration  count 

ISTATE  vector  for  E04KBF, 

workspace  for  VARIAB. 

If  covariance  matrices  are 
requested,  on  exit  IW0RK(2) 
contains  a  fault  indicator 
for  approximation  (i), 

IWORK( 3 )  contains  a  fault 
indicator  for  approximation  ( i 
and  IW0RK(4)  contains  the 
number  of  iterations  required 
to  compute  the  approximate 
Hessian  for  approximation  (i) 
Workspace  for  E04KBF ,  HESDIF 
Stores  MODE 

Stores  output  unit  number 
I  UNIT 

Stores  NPTS 

Stores  1  less  than  the 
address  of  PTS(l)  in  W 


PRECISION 

Double  precision  will  be  required  on  computers  with  32  bit  wordlength. 

The  precision  used  by  the  local  NAG  Library  implementation  should  be  ade¬ 
quate.  To  change  the  precision: 

-  change  all  REAL  declarations  to  DOUBLE  PRECISION; 

-  replace  constants  by  double  precision  versions,  constants  y  ,  y  , 

*10  typed  in  to  machine  accuracy; 

-  declare  NAG  Library  function  XQ2AAF  as  DOUBLE  PRECISION; 

change  the  precision  of  FORTRAN  library  functions,  i. e.  ,  ABS  to 

DABS,  ATAN2  to  DATAN2,  SIGN  to  DSIGN,  etc.  FLOAT(I)  can  be  replaced  by 

DBLE(FLOATCI)). 


If  extremely  large  observations  are  present  in  the  sample,  there 
may  be  a  loss  of  significant  figures  when  computing  sines  and  cosines 
in  SETECF  and  VMATRX'.  This  should  not  occur  when  real  data  is  used, 
but  can  be  a  problem  with  simulated  data  for  small  a. 


TIME 

Execution  times  depend  on  the  quality  of  initial  guesses  and 
properties  of  the  real  data  used,  and  vary  somewhat  throughout  the  para¬ 
meter  space.  As  a  rough  guide,  we  give  some  statistics  for  simulated 
data,  using  a  moderately  difficult  situation  with  a>l.  Tables  la  and  lb 
provide  approximate  running  times  for  Formulation  A,  q=40,  and  Formulation 
B,  q=20,  n=100, 200, 500, 1000, 2500.  Timing  starts  upon  entry  to  STABLE. 

Samples  from  S( 1.3,-. 5,3,15)  were  generated  using  the  method  of  Chambers, 
Mallows,  and  Stuck  (1976).  Initial  guesses  for  a,B,o,u  in  all  cases  were 
1505=^(1.01+2),  0,  ^(x  7,.-x  2^),  and  x  5,  the  sample  median,  respectively. 
Because  of  skewness,  the  median  is  not  a  good  estimator  of  u  in  this 
case.  Five  iterations  were  used.  Time  required  to  compute  asymptotic 
covariance  matrices  includes  approximations  (i)  and  (ii),  except  where 
noted.  Timings  are  for  a  double  precision  version  of  the  algorithm, 
compiled  by  the  IBM  FORTRAN  H  Extended  compiler,  and  run  on  an  IBM  370/3033. 

The  following  qualitative  points  are  clear  from  this  rather  restricted 
set  of  timings.  There  is  a  substantial  overhead,  which  may  crudely  be 
assumed  fixed,  associated  with  nonlinear  optimization,  although  E04KBF  solves 
the  optimization  subproblems  rapidly.  For  large  samples,  run  time  is 
dominated  by  evaluation  of  the  empirical  characteristic  function,  and 
thus  is  asymptotically  linear  in  n  for  a  fixed  number  of  iterations. 


Table  la 


Timings  for  Formulation  A, 
from  s(l. 3, -.5,3,15);  X=1  for 


n 

Iterations 

Estimation 
time  (sec) 

100 

5 

0.7 

200 

5 

1.0 

500 

5 

1.8 

1000 

5 

3.2 

2500 

5 

7.5 

q=40,  on  Simulated  Samples 


ns200  and  . 5 

for  n>200,  t=1 

Convergence 

Covariance 

criterion 

matrix  time 

5. 4( -4) 

0.1* 

2. 3( -2) 

0.1* 

1. 8( -4 ) 

0.8 

3. 3(  —  5 ) 

1.2 

4.  8( -6 ) 

2.5 

*Saraple  size  too  small  to  compute  approximation  (i),  only 
approximation  (ii)  computed. 


Table  lb 


Timings  for  Formulation  B,  q=20,  on  Simulated  Samples 


from  s( 1. 3,-. 5 

,3,15);  X=T=0 

Estimation 

Convergence 

Covariance 

n 

Iterations 

time  (sec) 

criterion 

matrix  time 

100 

5 

0.9 

1. 2( -2) 

0.3 

200 

5 

1.2 

3. 2( — 2 ) 

0.4 

500 

5 

1.6 

1. 8(  -4) 

0 . 5 

1000 

5 

2.3 

5. 7(  -5) 

0.7 

2500 

5 

4.4 

1.  5(  -4) 

1.4 
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Approximation  (ii)  of  the  asymptotic  covariance  matrix  is  quite  easy  to 
compute,  while  approximation  (i)  is  highly  time-consuming. 

For  fixed  k>l  with  the  k-sum  procedure,  one  iteration  generally 
suffices,  provided  estimates  from  the  nearest  value  of  k  are  used,  and 
the  estimates  don't  change  much.  For  mixtures  of  very  different  distri¬ 
butions,  or  if  the  exponent  a  is  near  unity,  more  are  required. 

ADDITIONAL  COMMENTS 

Although  output  need  not  be  produced,  we  recommend  calling  STABLE 
with  IUNIT>0,  so  the  user  will  have  a  record  of  how  estimation  progressed. 
The  following  information  will  then  be  printed  out: 

by  MONIT:  number  of  E04KBF  iterations  and  function  evaluations,  objective 
function  value,  norm  of  projected  gradient,  subproblem  solution, 
projected  gradient,  and  estimated  condition  number  of  projected 
Hessian; 

by  STABLE:  E04KBF  fault  indicator,  and  value  of  convergence  criterion; 

by  HESDIF( if  called) :  number  of  iterations  needed  to  compute  approximate 
Hessian,  and  step length  used. 

Use  of  STABLE  in  "batch  mode"  has  drawbacks.  For  instance,  most 
faults  arising  in  EQ4KBF  are  not  diagnosed  until  iteration  ceases.  In 
practice,  such  faults  may  likely  be  due  to  the  initial  a  being  on  the 
wrong  side  of  1.  Further,  when  a  is  small,  convergence  tolerances  are 
difficult  to  interpret,  and  the  user  may  prefer  direct  control  of  iteration. 
We  therefore  prefer  to  use  STABLE  interactively,  a  copy  of  the  output 
described  above  being  directed  to  the  terminal,  and  the  user  deciding 
after  each  iteration  whether  he  wishes  to  continue.  Required  modifications 
are  simple. 
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Faster  and/or  more  compact  codings  of  the  Algorithm  are  possible, 
for  instance,  if  8  is  known  to  be  zero,  if  only  Formulation  A  or  Formulation 
B  is  desired,  or  if  asymptotic  covariance  matrices  are  not  desired. 
Generality  is  achieved  at  a  price  in  efficiency. 
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CALL  FUNCT(0,  HPAR,  PAR,  TEMP,  PAR,  IW,  LIW,  W,  LW)  HESD  065 
PAR(IHD3)  =  PARJ  -  STEP  HESD  066 
CALL  FUHCT(0,  HPAR,  PAR,  TEMPI,  PAR,  IW,  LIW,  W,  LW)  HESD  067 
TEMP  =  TEMP  -  TEMPI  HESD  068 
PAR( I HD2 )  =  PARI  -  STEP  HESD  069 
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